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ABSTRACT 



Context. Discs in binaries have a complex behavior because of the perturbations of the companion star. Planetesimals growth and 
planet formation in binary-star systems both depend on the companion star parameters and on the properties of the circumstellar disc. 
An eccentric disc may significantly increase the impact velocity of planetesimals and therefore jeopardize the accumulation process. 
Aims. We model the evolution of discs in close binaries including the effects of self-gravity and adopting different prescriptions to 
model the disc's radiative properties. We focus on the dynamical properties and evolutionary tracks of the discs. 
Methods. We use the hydrodynamical code FARGO and we include in the energy equation heating and cooling effects. 
Results. Radiative discs have a lower disc eccentricity compared to locally isothermal discs with same temperature profile. Its 
averaged value is about 0.05, and it is almost independent of the eccentricity of the binary orbit, in contrast to locally isothermal 
disc models. As a consequence, we do not observe the formation of an internal elliptical low density region as in locally isothermal 
disc models. However, the disc eccentricity depends on the disc mass through the opacities. Akin to locally isothermal disc models, 
self-gravity forces the disc's longitude of pericenter to librate about a fixed orientation with respect to the binary apsidal line (n). 
Conclusions. The disc's radiative properties play an important role in the evolution of discs in binaries. A radiative disc has an overall 
shape and internal structure that are significantly different compared to a locally isothermal disc with same temperature profile. This 
is an important finding both for describing the evolutionary track of the disc during its progressive mass loss, and for planet formation 
since the internal structure of the disc is relevant for planetesimals growth in binary systems. The non-symmetrical distribution of 
mass in these discs causes large eccentricities for planetesimals that may affect their growth. 

Key words. Protoplanetary disks — Methods: numerical — Planets and satellites: formation 



1. Introduction 

Protoplanetary discs are known to exist in pre- main-sequence 
binary star systems through direct imaging ( Koerner et alj 
Il993t IStapelfeldt et al.lll99a iRodrfguez et al.lll998l) and spec- 
tral energy di stribut ions (lGhe 'z*etalTl l993l: iPrato et al ] 12003b 
iMonin et all |2007|) . iPetr-Gotzens et al.l (1201 Ol) have recently 



Of particular importance are the initial stages of planetesi- 
mal evolution when the mutual impact velocities must remain 
low to allow mass growth despite t he perturbations by the 



companion sta r . Prev i ous studies (e .g., Marzari & Scholll 
Thebault et all l2006l [2008L 120091: IPaardekooper et al] 



2000; 



2008 



Xie & Zhoul2009F Paardekooper & Leinhardtl2010l) have shown 



claimed that as much as 80% of all binary star systems in 
the ~ 1 Myr old Orion Nebula Clust er might have an ac- 
tive accretion disc. iTrilling et all J2007h showed that the inci- 
dence of debris discs in binary systems is not that different 
than that for single stars. These observations seem to suggest 
that the ground for planet formation in binary star systems is 
present even if, at subsequent stages, the gravitational perturba- 
tions by the companion star may neg at ively impact the growth 
proces s. iBonavita & Desideral d2007l) : iMugrauer & Neuhauserl 
(120091) have shown that the frequency of planets in binaries is 
not very different from that around single stars. However, this 
frequency critically depends on the binary separation and close 
binaries appear to be less favorable environments for planet for- 
mation. Only a few planets are presently known in tight systems 
(GI86, HD41004 and y Cephei) and it is expected that when the 
binary separation is in the range 20-100 AU the gravitational 
influence of the companion star must have influenced the planet 
formation process, affecting eit her the disc evolution or the plan- 
etesimal accumulation process dDesidera & Barbie"rill2007l) . 



that the combination of gas drag force and secular perturbations 
by the secondary star has strong effects on planetesimal orbits, 
leading to pericenter alignment and to an equilibrium distribu- 
tion for the eccentricity of small planetesimals. Since the magni- 
tude of the gas drag force depends on the planetesimals size, the 
variation of eccentricity and pericenter longitude with particle 
size may well inhibit the formation of larger bodies by exciting 
large mutual impact velocities. 

One question that needs to be addressed is to wha t extent 
the dis c' s perturbation du e to the companion star (e.g., ILubowl 
Il991allbl: iKlev et al.ll2008l) affects the distributions of eccentric- 
ity and longitude of pericenter for planetesimals of different size. 
Models in which the planetesimals orbital evol ution was com- 
puted along with the time evolu tion of the disc dKlev & Nelson! 
120081: IPaardekooper et al.ll2008l) have shown that the disc eccen- 
tricity may strongly influence the dynamics of small planetesi- 
mals by introducing radial gas drag forces and non-axisymmetric 
components in the gravity field of the disc. This additional per- 
turbing force may act either in favor or against planetesimal ac- 
cretion by affecting the relative impact velocities between the 
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bodies. An in depth analysis is needed to explore the strength 
of these perturbations. This first requires investigating the im- 
pact of the companion star on the disc's eccentricity. We fo- 
cus on relatively close binary star systems with large mass 
ratios and potentially large eccentricities which, according to 
iDuquennov & Mavorl(ll991l) . populate the peak of the frequency 
distribution in our neighborhood. Modelling the response of a 
protoplanetary disc to the gravitational perturbations by a sec- 
ondary (companion) star requires by necessity a numerical ap- 
proach, since it is difficult to predict analytically the disc shape 
and evolution with time. 

In a previous paper (iMarzari et al ] (120091) . hereafter Paper I) 
we used the numerical code FARGO to model the time evolution 
of a two-dimensional (2D) circumstellar disc in close binary star 
systems, including the effects of disc self-gravity. We focused on 
massive discs, with an initial mass equal to O.O4M , as well as a 
large binary's mass ratio (p = M s /M p = 0.4). This value is sta- 
tistic ally_the_2iiostfrec|^ binary systems observed 
so far (IDuquennov & Mavorlll991l) . We found that self-gravity 
significantly affects the increase in the disc eccentricity induced 
by the companion perturbations, and also influences the orienta- 
tion of the disc relative to the binary reference frame by causing 
libration instead of circulation. We also sampled different values 
of the binary's eccentricity, e^, ranging from to 0.6. The main 
findings of paper I were the following: 

- self-gravity plays a significant role in shaping the disc. The 
dynamical eccentricity ej (defined in Sect. 2.3) of the disc 
typically ranges from ~ 0.05 to ~ 0.15, depending on the 
binary eccentricity, e/,. It is smaller than in models without 
self gravity, 

- is inversely proportional to e/,, with the case of a circular 
binary {e/, = 0) being the most perturbing configuration, 

- the disc orientation ci> c / (defined in Sect. 2.3) librates around 
n, while it was circulating in the absence of self gravity, 

- an eccentric low-density region develops in the inner disc 
parts because of the large eccentricity and aligned pericen- 
ters of the gas streamlines there. 

The results of paper I have been obtained assuming a locally 
isothermal equation of state for the gas, wherein the initial ra- 
dial profile of the temperature remains constant in time, its value 
being set by the choice for the disc aspect ratio h - H/r, with 
H the pressure scale height. This approximation is well suited 
in disc regions where radiation cooling is efficient and the gas 
is optically thin. However, in particular in the initial stages of 
their evolution, discs may be optically thick and a more detailed 
treatment of the energy balance is required. In addition, when it 
passes at its pericenter, the secondary star triggers spiral shocks 
that may generate local strong shock and compressional heat- 
ing, which may violate the local isothermal approximation. The 
propagation of these shock waves may also be significantly al- 
tered in radiative discs. In this paper, we focus on how the disc 
eccentricity and orientation depend on the disc radiative prop- 
erties. Our approach is the following. We first examine how our 
previous results on locally isothermal discs depend on the choice 
for the (fixed) temperature profile of the disc. We then consider 
disc models with a radiative energy equation, for which we find 
that the averaged disc's eccentricity in its inner parts is smaller 
compared to locally isothermal disc models with similar temper- 
ature profile. We finally discuss the impact of our results in terms 
of planetesimals dynamics. 



2. Model description 

2.1. Code 

Two-dimensional hydrodynamical simulations were carried out 
with the ADSG versiorQ of the code FARGO. The code solves 
the hydrodynamical equations on a polar grid, and it uses an 
upwind trans port scheme alo ng with a harmonic, second-order 
slope limiter J van Leerll977h . The ADSG version of the FARGO 
code includes an adiabatic energy equatio n and a self-gravity 
modu le based on fast Fourier transforms dBaruteau & Massed 
I2008I) . Heating and cooling source terms have been implemented 
in the energy equation as described in Sect. I2.2I The specificity 
of the FARGO algorithm is to use a change of rotating frame on 
each ring of the grid, which increases the timestep significantly. 

Results of simulations are expressed in the following code 
units: the mass unit is the mass of the primary star, M p , which 
is taken equal to 1M . The length unit is set to 1 AU, and the 
orbital period at 1 AU is 2n times the code's time unit. 

2.2. Physical model and numerical setup 

We adopt a 2D disc model in which self-gravity and an energy 
equation are included (unless otherwise stated). The hydrody- 
namical equations are solved in a cylindrical coordinate system 
{r, ip] centered onto the primary star, with r e [0.5 AU - 15 AU] 
and tp e [0, 2n]. The grid used in our calculations has N r = 256 
radial zones and N s = 512 azimuthal zones, and a logarithmic 
spacing is used along the radial direction. The frame rotates 
with the Keplerian angular velocity at the binary's semi-major 
axis, and the indirect terms accounting for the acceleration of 
the primary due to the gravity of the secondary and of the disc 
are included. 

Binary parameters — Throughout this study, we adopt as 
standard model a binary system where the secondary star has a 
mass M s = 0.4M G . The binary is held on a fixed eccentric orbit 
with semi-major axis a\, = 30 AU and eccentricity e\, = 0.4, 
corresponding to an orbital period « 134 yr. 

Energy equation — Since the purpose of this work is to examine 
the impact of an energy equation on the disc's response to the 
periodic passages of a close companion, we carried simulations 
including either an energy equation, or a locally isothermal equa- 
tion such that the initial temperature profile remains constant in 
time. In all cases, the disc verifies the ideal gas law, 

P = KET/H, (1) 

where p and T denote the vertically-integrated pressure and tem- 
perature, respectively, E is the mass surface density, K is the ideal 
gas constant, and p is the mean molecular weight, taken equal to 
2.35. When included, the energy equation takes the form: 

JL + V • (ev) = -pV ■ v + g+ sc - Q- cml + AeV 2 log(p/l7) (2) 

where e = p/(y - 1) is the thermal energy density, y is the adi- 
abatic index, and v denotes the gas velocity. We take y = 1.4 
throughout this study. In Eq. (0, Q^ isc denotes the viscous heat- 
ing. We use both a constant shear kinematic viscosity, v = 10~ 5 
in code units, and a vo n Neumann-Richtmyer ar tificial bulk vis- 
cosity, as described in I Stone & Normanl dl992l) . where the co- 
efficient C2 is taken equal to 1.4 (C2 measures the number of 

1 See: http://fargo.in2p3.fr 
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zones over which the artificial viscosity spreads out shocks). The 
cooling source term in Eq. (0, Q moV is taken equal to 2os B T* s , 
where <tsb is the St efan-Boltzman n constant, and T e ff is the ef- 
fective temperature dHubenvlll990h . 



and 



eff 



r 4 /Teff, 



with effective optical depth 
3t V3 1 



(3) 



(4) 



The vertical optical depth, r, is approximated as r = kL/2, 
where for the R osseland mean opacity, k, the formulae in 
iBell & Lir] (1 19941) are adopted. Following iPaardekooper et afl 
J201 ll), we also model thermal diffusion as diffusion of the 
gas entropy, s, defined as s = 'Rly - 1) log(p/S y ). This 
corresponds to the last term in the right-hand side of Eq. (0, 
where A is a constant thermal diffusion coefficient. Throughout 
this study, we adopt A = 1CT 6 in code units. 

Initial conditions — The disc is initially axisymmetric and the 
angular frequency Q.(r) about the primary star is computed 
taking into account the radial acceleration due to the pressure 
gradient, the gravitational acceleration due to the primary 
star only and the self-gravity of the disk. The initial gas 
surface density, So, is set as IoO) = 2 (1 AU) x (r/l AU)~ 1/2 
with Eo(lAU) = 2.5 x 10~ 4 in code units. Assuming the 
mass of the primary star is 1M , this corresponds to setting 
the disc surface density at 1 AU to * 2.2 x 10 3 g cm -2 . 
Beyond 11 AU, £o(/) is smoothly reduced to a floor value, 
Sfloor = 10~ 9 = 4 x l(T 6 £o(l AU), using a Gaussian function 
with standard deviation equal to 0.4 AU. To prevent numerical 
instabilities caused by low density values or by steep density 
gradients near the grid's outer edge, the gas density in each 
grid cell is res et to Sfi nnr when ever it becomes smaller than 
this floor value dKlev et al.ll2.Q08l Paper I). Similarly, we adopt 
a floor value for the thermal energy density, efl 00 r = 10~ 18 . 
Our choice for this parameter is conservative, and we checked 
that with a larger efl oor , the modelled disc behaved similarly. 
The initial gas temperature, 7o, is taken proportional to r : 
T = 7o(lAU) x (r/l AU)" 1 , where the value for 7o(lAU), 
which is taken as a free parameter in our study, will be specified 
below. Writing the pressure scale height H as H — CsO K \ 
where cs denotes the sound speed and Ok the Keplerian 
angular frequency, the disc temperature can be conveniently 
related to the disc's aspect ratio, h — H/r. In our code units, 
where ft/fi = 1, this yields 7o(l AU) = 630Kx{/i(l AU)/0.05} 2 . 

Boundary conditions — By default, an outflow zero-gradient 
boundary condition is adopted at both the grid's inner and outer 
edges. The azimuthal velocity is set to its initial, axisymmetric 
value. No mass therefore flows back into the system, and the 
disk mass declines with time. The impact of the inner boundary 
condition on our results is examined in Sect. I4.3I 



2.3. Notations 

We will make use of the following quantities. We denote 
by and m& the disc's eccentricity an d perihelion longi - 



tude, respectively . They are defined a s in jKlev et all (l2008l) : 
iPierens & Nelson! d2007l );lMa 



larzari et 



ID d2009) 



M d 'x 



ff 



m& = M d 1 X J~J~ m(r, <p)%(r, <p)rdrdip 



(6) 



e(r, (fi)1.(r, <p)rdrdtp 



(5) 



where denotes the disc's mass, and e and vj are the eccen- 
tricity and pericenter longitude of each grid cell, respectively, 
assuming that the local position and velocity vectors uniquely 
define a 2-body Keplerian orbit. 



3. Locally isothermal disc models 

Before examining the impact of an energy equation on the disc's 
response to the tidal perturbations by the secondary star, we 
adopt in this section the simpler case for a locally isothermal 
equation of state, where the initial (axisymmetric) profile of the 
disc temperature remains fixed. This first step will help us an- 
alyze the more complex situation of a disc whose temperature 
evolves in time due to radiative cooling and various sources of 
heating, including that arising from the shock waves induced by 
the secondary. 



3. 1 . Disc's eccentricity and surface density profiles 

We carried out a series of simulations using a range of values for 
the disc's aspect ratio at 1 AU, h{ \ AU). This comes to varying 
the disc's temperature at the same location. Other disc and binary 
parameters are as described in Sect. 12.21 We consider two differ- 
ent values for the secondary-to-primary mass ratio: q = 0.4 (our 
fiducial value), and q — 0.1. Results of simulations with q = 0.4 
are shown in the upper panels in Fig. [TJ and those with q = 0.1 
in the lower panels. Azimuthally- and time-averaged profiles of 
the disc's eccentricity and surface density are displayed in the 
left and right columns of Fig. [TJ respectively. Profiles are dis- 
played after ~ 60 orbits of the secondary, and time averaging is 
done over 5 orbits. By checking the time evolution of the disc's 
eccentricity profile, we find that a steady state is reached after 
typically 20 binary revolutions with q = 0.4 for all values of h, 
and after 40 revolutions for q = 0.1. The value of h(l AU) is 
indicated in each panel (and simply denoted by h). 

We are primarily interested in the disc's averaged eccen- 
tricity in its inner parts, below the truncation radius located at 
r ~ 5 - 6 AU for q = 0.4, where planet forma tion is more 
likely to occur. As shown in iMarzari eta D d2009h . the trunca- 
tion radius of the disc closely matches the critical limit for or- 
bital stability due to the secula r perturbations of the companion 
star dHolman & Wiegertlll999h . From Fig. Q] it is clear that the 
disc eccentricity increases from r ~ 4 AU downwards, peaks at 
r < 1 AU, then decreases towards the location of the grid's inner 
edge, where the boundary condition imposes zero eccentricity. 
Interestingly, we see that the averaged peak eccentricity of the 
disc (reached at r < 1 AU) increases with /; up to /; = 0.05 and 
decreases beyond this value. The same behavior is obtained with 
both values of q. This behavior may be interpreted as follows. 
The disc's density perturbation due to the secondary decreases 
with increasing disc aspect ratio (h or, equivalently, increas- 
ing sound speed). In the limit of large aspect ratios, the disc's 
eccentricity thus decreases with increasing h. As h decreases, 
the disc's perturbed density increases, but shock waves become 
more tightly wound, and have to travel a longer distance before 
reaching the disc's inner parts. Being then more prone to viscous 
damping, shock waves become less and less efficient at deposit- 
ing angular momentum in the disc's inner parts, where the ec- 
centricity remains small. This explains why when decreasing 



3 



F. Marzari et al.: Radiative discs in close binaries 



the averaged eccentricity of the disc's inner parts first reaches a 
maximum and then decreases. 

We can see in Fig.[T]that the radial dependence of the den- 
sity and eccentricity profiles reasonably match each other, as 
expected. The general trend is that the larger the disc's eccen- 
tricity profile, the smaller the surface density profile. Note that 
this is more visible for q = 0.4, where the density's perturba- 
tion by the secondary takes larger values than for q = 0.1. The 
significant decrease in the disc's density profile in the disc inner 
parts (typically below ~ 2 AU) takes the form of an elliptic inner 
hole, which has been observed in a number of previous numeri- 
cal studies (e.g., lKlevetaill2008L Paper I). 

Our results concerning the dependence of the di sc's eccen- 
tricity and density profiles with h differ from those of iKlev et al.l 
d2008l) . but this depends on intrinsic differences in the models. 
In particular, we consider more massive discs and we include the 
effects of self-gravity which proved to be efficient in affecting the 
disc evolution. 

3.2. Fourier analysis of the density distribution 

To provide some insight into the dependence of the disc's eccen- 
tricity and density profiles with varying temperature, we exam- 
ine in this paragraph the Fourier components of the disc's sur- 
face density. Fig. [3] displays the instantaneous amplitude of the 
Fourier mode coefficients with azimuthal wavenumber 1 < m < 
5. Results are shown at 2550 yr, that is after « 20 orbits of the 
secondary, when the disc is truncated at about 7-8 AU. We com- 
pare the profiles obtained for previous series of locally isother- 
mal disc models for q = 0.4, with /; = 0.02 (left panel) and 
h = 0.06 (right panel). The secondary star is half-way between 
the pericenter and apocenter at this particular point in time. 

From Fig. |3] it is clear that the amplitude of Fourier mode co- 
efficients decreases with increasing m, and that the m = 1 mode 
prevails. The run with h = 0.06 displays larger mode amplitudes 
throughout the disc's inner parts, which accounts for the larger 
disc eccentricity obtained in this case. Also, the large amplitude 
of the m = 1 mode below r m 1 for h = 0.06 is directly associ- 
ated to the presence of an elliptic inner hole, whose presence has 
already been pointed out in the upper panels in Fig. Q] 




50 - 

I i i i i i i i 
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R (AU) 

Fig. 2. Azimuthally-averaged profile of the disc's pericenter in 
the locally isothermal disc model with h = 0.05 and q = 0.4. It 
is computed after 100 orbits of the secondary star. 



4. Radiative discs models 

In this section, we describe the results of our hydrodynamical 
simulations that include an energy equation. Our aim is to assess 
the impact of the energy equation on the disc's averaged eccen- 
tricity. 

4. 1 . Disc eccentricity 



'o 
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Fig. 5. Time evolution of the averaged disc eccentricity ej (up- 
per panel) and pericenter longitude Wd (lower panel) for the lo- 
cally isothermal and radiative runs of Sect. I4.ll 

As a first step, we compare a radiative and a locally isother- 
mal model with same initial aspect ratio, h = 0.05. All other disc 
parameters are as described in Sect. [2] Fig.|4]shows contours of 
the disc's surface density at 90 binary orbits obtained with two 
values of the binary's eccentricity: et, = 0.4 (upper panels) and 
e/, = (lower panels). As already pointed out in Sect. [3] eccen- 
tric streamlines (m = 1 density mode) have different values of 
the pericenter depending on the radial distance, and, as a conse- 
quence, they combine into a pattern of spiral structure. However, 
the discs computed with the radiative model (left-handed plots 
in Fig. |4|i appear smoother and more symmetric than the corre- 
sponding locally isothermal discs, and, in particular, they do not 
feature any elliptic hole near the inner edge, independently of 
et,. The spiral density waves are less strong in the radiative case 
and this is possibly related to the absence of a low density region 
close to the star. The disc eccentricity e,; and perihelion longitude 
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02468 10 02468 10 

R (AU) R (AU) 

Fig. 1. Azimuthally- and time-averaged profiles of the disc's eccentricity (left column) and surface density (£, right column) obtained 
with the series of locally isothermal disc models described in Sect. [3] In the density profile plots the initial, unperturbed profile is 
shown as reference. Results are shown at 60 orbits of the secondary, and time-averaging is done over 5 orbits. Several values of the 
disc's aspect ratio at 1 AU are considered and two secondary-to-primary mass ratios: q = 0.4 (upper panels) and q = 0.1 (lower 
panels). 



are shown as a function of time in Fig. [5] and they confirm that 
the radiative case has, on average, a lower eccentricity even if it 
takes more time to reach a steady state. The perihelion libration, 
observed also for the radiative case, strongly suggests that this 
behaviour is solely due to the disc self-gravity and that it does 
not depend on the energy equation. By inspecting the radial peri- 
center profile, even in the radiative case the azimuthally averaged 
pericenter of the gas streamlines changes with radial distance (as 
is illustrated in Fig|2]for a locally isothermal model). This vari- 
ation, in addition to causing spiral waves, leads to an asymmet- 
ric distribution of mass and then to a non-symmetric disc gravity 
field. This is critical forplanetesimals embedded in the disc since 
the non-homogeneous disc forces non-radial components on the 
gravity field felt by planetesimals significantly perturbing their 
orbits. These perturbations are comparable in magnitude to the 
secular effects of the companion star but are irregular and may 
then cause larger changes in the planetesimals orbital elements. 
They are indeed an indirect effect of the companion gravity but 
for planetesimals they represent an independent source of per- 
turbation. 

To get further insight into the different densities and eccen- 
tricities with and without an energy equation, we carried out 
two additional locally isothermal disc models. In the first model, 
the initial temperature is set to the time-averaged temperature 



profile of the above radiative run. In the second model, the ini- 
tial temperature is chosen to give the same sound speed profile 
as the time-averaged one in the radiative run. In Sect. 13.11 we 
pointed out that the disc eccentricity in locally isothermal mod- 
els strongly depends on the aspect ratio h and, as a consequence, 
on the temperature profile (see Fig. Q}. In Fig. [6] bottom plot, 
we compare the temperature profile of the locally isothermal 
model with h = 0.05 to the time-averaged temperature profile 
of the radiative model. The temperature in the radiative run is 
significantly larger, implying that we should compare the radia- 
tive model to a locally isothermal model with h ~ 0.07. In the 
top panel of Fig. [6] we display the eccentricity and density pro- 
files in all four models: (i) the radiative model, (ii) the locally 
isothermal model with same initial aspect ratio as in the radiative 
model, (iii) the locally isothermal model with same temperature 
profile as in the radiative case, and (iv) the locally isothermal 
model with same sound speed profile as in the radiative one. The 
disc eccentricity in models (iii) and (iv) is smaller than in model 

(ii) , as expected from the results shown in Fig.Q](top-left panel). 
Consequently, the disc surface density remains larger in models 

(iii) and (iv) compared to in model (ii). Still, models (iii) and 

(iv) do not match the low eccentricity of the radiative case and 
its smooth density distribution close to the star. The strength of 
the m — 1 Fourier mode in the radiative run is smaller in the 
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r(AU) r(AU) 

Fig. 3. Instantaneous Fourier components of the disc's surface density for the locally isothermal disc models of Sect.[3]with h = 0.02 
(left panel) and h = 0.06 (right panel), respectively. Results are shown at t — 2550 yr (that is, after about 20 orbits of the secondary), 
and are obtained with q = 0.4. The azimuthal wavenumber for each Fourier component is indicated in the top-right corner of each 
panel. 



disc's inner parts, and this difference depends on the form of the 
energy equation. In Fig. [7] we display density contours obtained 
with the radiative run (left panel), and with the locally isothermal 
model with same temperature profile (right panel). In the radia- 
tive model, the wave perturbations are less strong and the waves 
appear to be more damped. Being weaker in the inner disc parts 
in the radiative case, waves should be less efficient at deposit- 
ing angular momentum there. This would explain why radiative 
discs are less eccentric in the regions close to the central star, 
and why the surface density is more homogeneous. 

4.2. Radiative damping of waves 

The main question arising from previous results is what causes 
a stronger damping of density waves in the radiative model than 
in the locally isothermal one, while the same (time-averaged) 
sound speed profile is adopted in both models. In Fig. [6] the 
radiative model shows indeed a significantly lower eccentricity 
profile. Assuming the disc eccentricity is related to the strength 
of the density waves induced by the binary gravitational per- 
turbations, we have shown in Fig. [7] that spiral waves are more 
damped in the radiative case. 

There are two potential sources of wave damping that may 
help understand why the disc e ccentricity remains lower in the 
radiative case: sho ck damping (|Goodman & Rafikovll2001l) and 
radiative damping (ICassen & Woolumll 19961) . To explore the ef- 
ficiency of the shock damping mechanism, we examined the 
vortensity distribution in the disc, since it experiences a jump 
at shocks, the magnitude of which depends on the strength of 
the shock. We did not observe significant differences in the disc 
vortensity distribution between the two models and thus we be- 
lieve that shock damping does not significantly contribute to the 
smoother behaviour of the radiative disc. 

The energy los s by radiation is an additi onal possible source 
of wave damping dCassen & Woolum|[T996h . Wave propagation 
through adiabatic compressions and expansions may be damped 
by radiative losses which, in our model, are controlled by the 
cooling term Q~ r To test this hypothesis, we restarted our stan- 
dard radiative simulation (e/, = 0.4, M s = O.4M and «b = 30 
AU) after 165 binary revolutions adopting different two cooling 
prescriptions. In a first run, we limit the cooling time throughout 



the disc to be no less than that at 6 AU from the star (that is, 
about 150 yr). To keep the temperature profile as close as pos- 
sible to that of the standard radiative model, and thus to prevent 
discrepancies related to different temperature profiles, we limit 
the viscous heating timescale by a similar amount. In a second 
run, we increase the cooling rate throughout the entire disc by 
a factor of 10, and the viscous heating rate is increased accord- 
ingly. Both restart simulations were run for 50 additional binary 
orbits, over which the disc's temperature profile does not evolve 
significantly, except within 1 AU from the central star. 

Fig.[8]compares the disc time-averaged eccentricity and tem- 
perature profiles of the standard radiative model with those of the 
additional models with longer and shorter cooling times. The in- 
ner disc's eccentricity is significantly increased with longer cool- 
ing timescale, exceeding ~ 0.1 almost uniformly from R = 2 
AU to R = 6 AU. This value is in good agreement with that of 
the locally isothermal run with same sound speed profile (see 
Fig . [6j> . However, close to the inner edge the disc eccentricity is 
still small but this is possibly due to the fact that we cannot main- 
tain constant the temperature profile within 1 AU of the central 
star in the model with reduced cooling. In the restart simulation 
with a cooling rate 10 times larger, the disc eccentricity is about 
half that of the standard radiative model, while the correspond- 
ing temperature profiles hardly differ. 

The results shown in Fig.[8]suggest that radiative damping is 
responsible for the limited growth in the disc eccentricity com- 
pared to to similar locally isothermal disc models. This mecha- 
nism, like self-gravity, helps maintaining the disc eccentricity to 
a mild value. 

4.3. Dependence on the boundary conditions 

We examine in this paragraph the dependence of our results on 
the choice for the boundary condition at the grid's inner edge. 
For this purpose, we compare the results of simulations using 
our fiducial boundary condition (zero-gradient outflow bound- 
ary condition, see Sect. 14. U with those of two additional simu- 
lations: 

- One using a visco us outflow boundary condition 
dKlev & Nelson! 120081) . It is very similar to our zero- 
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plots are for = 0.4 and the lower plots for et, = 0.0. The same initial aspect ratio (h = 0.05) is used in these simulations. 



gradient outflow boundary condition, except that the 
radial velocity at the inner boundary is set to the local 
(azimuthally-averaged) viscous inflow velocity of a disc 
in equilibrium with locally uniform surface density profile 
(-3v/2R m i n ). The azimuthal velocity is also set to its initial 
axisymmetric value, 

Another simulation also using our standard outflow bound- 
ary condition, but where the azimuthal velocity at the inner 
boundary is extrapolated from that in the first active ring with 
a r 1 ' 2 law. In contrast to previous boundary conditions, the 



disc at the grid's inner edge is no longer forced to remain 
circular. 

The results of this comparison for locally isothermal runs 
with same temperature profile as in the radiative model are 
shown in Fig. [9] The profiles are in good agreement for all three 
different boundary conditions. Note that the surface density in 
the inner disc for the viscous boundary condition is larger than 
with the two other boundaries. This is expected, since the im- 
posed viscous drift velocity is smaller that the radial velocity set 
by the propagation of the spiral waves induced by the secondary. 
A similar good agreement is obtained with a radiative model, as 
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illustrated in Fig.[l0] These results confirm the robustness of our 
results against the choice for the inner boundary condition. 



4.4. Dependence on the binary eccentricity 

In paper I, we showed that the disc eccentricity decreases with 
increasing binary eccentricity, the circular case (e;, = 0) being 
the most perturbing configuration. We interpreted this result as 
being due to the larger size of the disc for lower values of ei„ 
and to the consequent larger number of resonant perturbations 
which may affect it. The eccentricity of radiative discs, however, 
seems to be rather insensitive to <?/,. The outcome of the simu- 
lations at different e/, is shown in the left plot of Fig.Qj] where 
we compare the values of obtained for locally isothermal and 
radiative discs. In contrast to locally isothermal discs, the aver- 
aged of radiative discs is almost constant around 0.05, and it 
does not show a significant dependence with <?/,. However, the 
azimuthally-averaged profile of the disc's eccentricity over the 
radiative disc does depend on e/,. In Fig.[12j the radial profile of 
the disc's eccentricity is shown for different values of et, rang- 



ing from to 0.6. We note that the eccentricity profiles signif- 
icantly differ even if the position of the star with respect to the 
initial reference frame is the same. All radiative runs show an 
increasing eccentricity profile towards large radii, whereas lo- 
cally isothermal runs feature a peak in the eccentricity in the 
inner disc parts for most values of h. It should be pointed out 
that the above comparison is done by adopting the same initial 
aspect ratio for radiative and locally isothermal disc models. A 
different approach would be to compare disc models with same 
temperature or sound speed profiles, as is done in Figs. [7] and [6] 
But, when comparing the temperature profiles of the radiative 
discs, we find that they correspond to locally isothermal discs 
with h ~ 0.06 - 0.07, which anyway feature a larger eccentricity 
in the disc inner parts. 

In the right panel of Fig. QT|we also compare the values of 
the outer semi-major axis a c t of the ellipse best fitting the outer 
edge of the disc (at a density level of S = 10~ 6 , code units). 
The radiative discs extend slightly farther out compared to the 
corresponding locally isothermal discs, except for our fiducial 
binary's eccentricity (e^, = 0.4), where the disc's size is approx- 
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Fig. 11. Disc eccentricity (left plot) for different values of the binary eccentricity e\,. The average of ed is displayed after 70 binary 
revolutions, and time-averaging is done over 20 revolutions. For com parison, we also depict the corresponding values of e,/ for 
isothermal discs with same binary configurations dMarzari et al. 2009). In the right plot, we show the disc's size for different values 
of eb, more specifically the semi-major axis of the ellipse that best fits the disc density at a level X = 10~ 6 (code units). 



imately equal to its isothermal counterpart. The mass loss rate 
for radiative discs, for different values of <?£,, is shown in Fig. [13] 
It is computed with a linear fit of the time evolution of the disc 
mass when a steady state is reached, so it is independent from the 
initial disc truncation process. It shows an almost linear depen- 
dence with eb and it is related to the strong perturbation effects 
on the disc as the secondary star passes at pericenter. The inspec- 
tion of the amount of mass stripped vs. time reveals that almost 
all the mass loss occurs during and after the pericenter passage of 
the companion while a lower amount is constantly lost through 
the inner border due to the disc eccentricity and viscosity. The 
absolute value of the mass loss is large, and it predicts a reduc- 
tion in the disc mass by a factor of two in about 3 x 10 4 yrs for 
our standard case {eb = 0.4). However, this value depends on 
the disc mass, and a simulation with the same binary parameters 
but an initial disc mass 10 times smaller gives a mass loss rate 
* 3.7 x 10~ 9 M /yr. 

4.5. Dependence on the disc mass 

In Paper I, we found that for locally isothermal discs, a decrease 
in the disc mass had no significant impact on the disc eccen- 
tricity. Below we show that radiative discs behave differently. In 
Fig. [T4J we compare the disc eccentricity in the nominal case, 
where the initial density at 1 AU is that of the MMSN, to that of 
a disc initially 2 and 10 times less massive. The disc eccentric- 
ity ea in a steady state is much larger for the less massive discs. 
This is further illustrated by the contours of the disc density for 
Eo = O.IEmmsn in Fig. [15] The disc appears smaller and very 
eccentric, particularly in its outer regions. 

It is difficult to interpret this outcome on the basis of the 
isothermal simulations. Less massive radiative discs, like that 
shown in Fig. [15] have lower temperature profiles. This is pre- 
sumably due to a faster cooling rate. According to our model, a 
lower surface density of the disc implies a smaller optical thick- 
ness that leads to a shorter cooling timescale. As a consequence, 
our standard model with Eo = 1 Smmsn has a stationary tempera- 
ture profile which is equivalent to that of a locally isothermal run 
with h ~ 7 - 8% while less massive discs have a lower temper- 
ature profile in a steady-state with a temperature typically 3-4 
times less high. This translates into equivalent locally isother- 
mal models with h down to 3 - 4%. By inspecting Fig. Q] (top 
left plot) we notice that in this range of aspect ratios, the disc's 
eccentricity actually changes quickly with This might explain 



why less massive radiative discs have different disk eccentricity. 
Even self-gravity will play a minor role for less massive discs 
and in Paper I we showed that indeed self-gravity is effective 
in reducing the disc eccentricity. This, however, is not the full 
story and the different energy equation also plays a significant 
role. The disc eccentricity radial profile of radiative discs, also 
the less massive ones (see Fig. fT4~t . does not peak close to the 
star like for the isothermal discs (Fig. []]). The value of ed for 
radiative discs grow s for larger values of R as predicted by the 
analytical theory of iPaardekooper et al.l d2008l) . This behaviour 
is further illustrated by examining the Fourier components of the 
normalized surface density. Fig. Q~6] compares the Fourier com- 
ponents for both our standard disc with So = Emmsn, an d me 
disc model with So = O.ISmmsn- The m = 1 and m = 2 com- 
ponents are much stronger for the less massive disc, and they 
account for the overall higher disc eccentricity. In the standard 
case, the high disc density is able to damp efficiently the two 
Fourier components of the binary perturbations as they move to- 
ward the disc inner parts. On the other hand, these components 
propagate further in for less massive discs, leading to a larger 
value for e,/. However, also in the less massive disc they do not 
reach the inner part of the disc, like in isothermal discs, prevent- 
ing the formation of an eccentric low-density region close to the 
central star. 

According to the results presented so far, isothermal and ra- 
diative discs behave differently. Massive isothermal discs in low 
eccentricity binaries are expected to be eccentric for reasonable 
values of /; while radiative discs always have low eccentricity. 
For isothermal discs there is a weak dependence of the disc 
eccentricity on the disc mass (see Paper I) while, for radiative 
discs, this dependence is strong with less massive discs being 
more eccentric. Hot isothermal discs develop an inner eccentric 
hole while radiative discs are smooth also at the inner edge. 

5. Why is disc eccentricity so important? 
Planetesimals! 

The shape and profile of the gas disc may have a crucial role 
in the early stages of planetary formation, when kilometer-sized 
planetesimals are colliding with each other. Even when neglect- 
ing gas disc gravity, several studies have shown that the coupling 
between secular perturbations and gas drag, which increases im- 
pact velocities between non-equal s ized bodies and can lead to 
accretion-hostile environments (e.g. lThebau lt et al. 2 0061 [2008L 
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Fig. 12. Azimuthally-averaged eccentricity profiles for radiative 
discs with different binary eccentricities, displayed at 100 binary 
orbits. Instantaneous profiles are depicted to better highlight the 
dependence of with eb- 



Fig. 13. Mass loss rate of radiative discs (in M e /yr) for different 
values of obtained through a linear fit of the time evolution of 
the disc mass after the initial fast truncation. 
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Fig. 14. Azimuthally- and time-averaged profiles of the disc's 
eccentricity for three different initial surface densities: So = 
Smmsn (our fiducial value), So = 0.52 M msn, an d So = 
O.ISmmsn- 



l2009h . leads to further increasing the relati ve velocity between 
plane t esimals if the g as disc gets eccentric dPaardekooper et all 
120081 : IXie et al.ll2Q10h . Disc gravity might further increase im- 
pact velocities by increasing the disc eccentricity and inducing 
additional dynamical perturbations on planetesimals because of 
non-isotropic distribution of the mass within the gas disc. The 
asymmetric distribution of mass within a massive disc can in- 
deed significantly perturb the orbit of planetesimals, causing 
large eccentricities and unphased orbits, which may possibly halt 
the accretion process. 

To test this hypothesis we integrated the trajectories of plan- 
etesimals orbiting within 6 AU from the primary star, and 
we computed their orbital evolution in 2D under the action 
of (i) stellar gravity, (ii) gas drag, and (iii) gas disc grav- 
ity. This procedure has been implemented in previous papers, 
but with different assumptions like isother mal disc models, no 
self-gravity and different bin ary parameters dPaardekooper et all 
12008b iKlev & Nelson! 120081 We first present in Fig. [17] the re- 
sults of a test run without gas drag, for which we see that the 
planetesimals eccentricity grows to large values regardless of 
their initial location in the disc. In addition, the drift towards 
the inner region of the disc, which is typical of planetesimals 
around single stars, may be halted and even reversed. This drag- 
free case could probably describe the evolution of large plan- 
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Fig. 16. m = 1 and m = 2 Fourier components of the normal- 
ized surface density for a radiative disc with So = ^mmsn, and a 
radiative disc initially ten times less massive (Zq = O.ISmmsn)- 



etesimals, 100 km or bigger, which are not significantly affected 
by gas drag. For the disc and star parameters, we consider our 
standard case. 

We then present the results of a full simulation including 
gas drag, whose expression is given by the usual formulae Fa 
K | v le i I v re i, where v re i is the relative velocity vector of the plan- 
etesimal with respect to the gas ; and the drag parameter K is 
equal to p g Cd/(8p p is) dKarv et al.l ll993j). where s is the radius of 
a given planetesimal, p p \ its mass density, p g the gas density of 
the protoplanetary disc, and Cd a dimensionless drag coefficient 
related to the planetesimals shape (~ 0.4 for spherical bodies). 
We explore 3 values for K, corresponding to 10, 20 and 50 km- 
sized planetesimals in a MMSN disc, and we consider 2 initial 
locations for the planetesimals: 1.5 AU and 3.5 AU from the pri- 
mary. For the 1.5 AU case, the initial strong perturbations of the 



eccentric disc are partly damped by gas drag but the steady state 
eccentricities are much larger than the forced eccentricity of the 
companion star (see Fig. [18] upper box). Moreover, these steady 
state eccentricities do strongly vary with K, i.e. with planetesi- 
mal sizes: they are almost twice as large for 50 km-sized objects 
as for 10 km-sized ones. As for the pericenter longitudes, they 
rapidly converge to steady state values that only marginally vary 
with planetesimal size. The evolution of the 3.5 AU case is dif- 
ferent. The eccentricities decrease towards an equilibrium value 
predicted by the balancing between the forced component of the 
companion star, gas drag damping and disc gravity (see Fig. [18] 
lower plots) The steady state eccentricities are thus lower than in 
the 1 .5 AU case, despite being in a region where secondary per- 
turbations are st ronger. This result is in sha rp contrast with what 
was obtained bv lPaardekooper et al.l (1200 8) for the gas drag-only 
case, where eccentricities steadily increase with increasing semi- 
major axis (see Fig. 10b of that paper). It clearly illustrates the 
fact that the gas disc gravity, which is stronger in the inner and 
denser part of the system, is the dominant mechanism control- 
ling the planetesimals' dynamical evolution. The residual short- 
term variations of the eccentricity, due to the companion's per- 
turbations, are in contrast much larger in the 3.5 AU case than in 
the 1.5 AU one. The pericenter alignment is maintained at this 
distance even if it is less collimated, as expected since the gas 
density is lower. 

Assessing the consequences of these dynamical behaviors 
on the accretion process would require estimating the distribu- 
tion of impact velocities, v co n, among the planetesimals popula- 
tion. It is here not possible to directly derive v co \\ from the val- 
ues of the eccentricity because of mutual orbital phasing. The 
simple v co n oc (e) relation is in this case no longer valid and 
should be replaced by v co n oc ap (e), where the factor ap < 1 
accounts for the phasing. Unfortunately, proper velocity esti- 
mates are difficult to derive with the limited number of plan- 
etesimals (50) considered here. Simulations wit h at least a few 
thous a nds test particles would be necessary (e . g.lThebault etail 
120061: iPaardekooper et al.l 120081: IXie & Zhoul l2009h . but such 
numbers are beyond the current computing capacities for a full 
model including the gas disc's gravity. However, the preliminary 
results displayed in Figs. [17] and [18] seem to indicate that en- 
counter velocities are probably further increased by the action 
of disc gravity, and this for 2 reasons: (i) Even if a significant 
orbital phasing is observed (right hand plots of Fig. [T81 it is 
nevertheless not perfect, especially in the 3.5 AU case. Even in 
the 1.5 AU case it is size-dependent, which increases the fac- 
tor ap in any realistic planetesimal population with a size dis- 
tribution, (ii) The values of the steady-state eccentricities are 
higher than in the gas-drag only case, especially closer to the 
star where velocities are relatively low in the gas-drag only runs 
(c ompare for example the up per-left graph of Fig. [18] to Fig. 10 
of IPaardekooper et al.l (120081) ). These steady-state eccentricities 
are in addition strongly size-dependent. These two effects sep- 
arately increase each of the two components, ap and (e), con- 
trolling the value of the encounter velocities. The extent of this 
velocity increase, and thus its concrete effect on the accretion 
process, cannot be quantitatively estimated here. It will be the 
purpose of a forthcoming study specifically addressing this is- 
sue. We are however confident that the general tendency is for 
disc gravity to act against planetesimal accretion. 

6. Concluding remarks 

We have shown that the evolution of a circumstellar disc in 
close binary-star systems strongly depends on the disc's mass 
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Fig. 17. Planetesimal orbital evolution without the gas drag 
force. Large values of eccentricity are excited by the non- 
homogeneous distribution of mass within the disc. The highest 
eccentricity value is observed closer to the star where the com- 
panion gravitational perturbations are weaker. 



and thermodynamical properties. In locally isothermal disc mod- 
els, where the temperature profile remains constant with time, 
the averaged eccentricity in the disc inner parts changes dra- 
matically with varying the disc temperatures (or, equivalently, 
the disc aspect ratios h = H/r), as illustrated in Fig. [TJ In our 
model, the disc eccentricity typically peaks at « 0.4 for h as 0.05. 
Such large eccentricities result in the formation of an internal 
elliptic low-density region in the disc. Radiative discs, on the 
other hand, have a smoother density profile, and their eccentric- 



ity takes smaller values than in locally isothermal models with 
same temperature profile. Radiative damping of the waves in- 
duced by the secondary companion contributes to keep the disc 
eccentricity to a fairly small value, which in our model typically 
amounts to ~ 0.05 in the disc inner parts. 

In both locally isothermal and radiative discs, self-gravity 
causes the libration of the disc orientation at an angle n with 
respect to the apsidal line of the binary orbit. However, the libra- 
tion is not coherent within the disc, and the orientation of the gas 
streamlines changes with radial distance. 

The averaged disc eccentricity in radiative discs is almost 
insensitive to the binary's eccentricity, although it should slowly 
increase with time as the disc mass decreases due to its viscous 
evolution. A disc that is 10 times less massive than our nominal 
disc (that is, having So = O.IZmmsn) has a large eccentricity 
(about 0.3 for e;, = 0.4). This behavior is not observed in locally 
isothermal discs. 

Based on these different behaviors, we may envision an evo- 
lutionary track for discs in binaries. In the earlier stages, the disc 
is massive and probably optically thick so it is well described by 
radiative models. Its eccentricity should therefore remain small, 
around 0.05, and its radial density profile be smooth. Later, be- 
cause of the viscous evolution, the disc progressively loses mass. 
This occurs at a fast rate for binaries with large values of e\, (see 
Fig.QT). If the disc is still optically thick, and the radiative model 
is appropriate, its eccentricity is expected to grow because of the 
dependence of e,/ with the disc mass (see Fig. [T4l . If, finally, 
the disc becomes optically thin, a locally isothermal equation of 
state may then be more appropriate. In this case, t he disc eccen- 
tricity will change and depend on e\, as shown in iMarzari et alj 
(120091) . At the same time, hotter discs will develop an internal 
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elliptic hole, a significant decrease in the gas density at the inner 
edge of the disc. 

The non-symmetric distribution of mass caused by the dif- 
ferential libration of the disc orientation and the disc overall ec- 
centric shape causes significant perturbations to planetesimals. 
Their eccentricity is driven to values significantly larger than 
those caused only by the secular perturbations of the binary com- 
panion. Even when gas drag is included, the perturbations by the 
eccentric disc entail large planetesimals eccentricities, in partic- 
ular in the inner zones of the disc. It is necessary to evaluate the 
impact of this eccentricity on the accumulation process of plan- 
etesimals. 
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